1a - Apply PCA, NMF, ICA and MDS, UMAP, and tSNE to this dataset. Compare and contrast the results using these methods.
1b - Relate the dimension reduction results with the clinical data. Is any clinical information reflected in the lower dimensional spaces?
1c - Overall, which dimension reduction method do you recommend for this data set and why?
2a - Apply various clustering algorithms such as K-means (explore different K), hierarchical clustering (explore different linkages), NMF, and biclustering. Compare the clustering results using these methods.
2b - Relate the clustering results with the clinical data. Can the clustering algorithm recover some of the clinical information such as cancer subtypes?
2c - Use Consensus Clustering to help Validate Clustering Results
2c - Overall, which clustering method(s) do you recommend for this data set and why?
3a - Identify important genes to differetiate ER postive and negative, PR postive and negative, HER2 postive and negative, metastasis status.
3b - Try different procedures to adjust for multiple comparisons.
3c - Examine the lists of genes identified using different methods for each clinical response. Which method is best? Why?
5a - Use graphical models to explore interactions among genes. Are any of the well-connected genes related to patterns previously identified?
6a - Visualize this data using multiple approaches.
6b - Prepare the “best” visual summary of this data.
7a - What is the most interesting finding?
7b - Is this finding consistent and stable?
7c - Prepare a visual summary that best illustrates this interesting finding.
R scripts to help out with the BRCA case study Lab Don’t peek at this if you want to practice coding on your own!!
Load Data
load("UnsupL_SISBID_2020.Rdata")
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 3.6.2
library(ConsensusClusterPlus)
library(kknn)
library(GGally)
## Warning: package 'GGally' was built under R version 3.6.2
## Registered S3 method overwritten by 'GGally':
## method from
## +.gg ggplot2
library(umap)
## Warning: package 'umap' was built under R version 3.6.2
library(Rtsne)
library(igraph)
##
## Attaching package: 'igraph'
## The following objects are masked from 'package:stats':
##
## decompose, spectrum
## The following object is masked from 'package:base':
##
## union
library(XMRF)
library(spacejam)
## Loading required package: splines
## Loading required package: Matrix
Explore Data
dim(gdat)
## [1] 445 353
dim(cdat)
## [1] 445 6
table(cdat$Subtype)
##
## Basal-like HER2-enriched Luminal A Luminal B Normal-like
## 79 53 200 106 7
table(cdat$ER)
##
## Indeterminate Negative
## 2 100
## Not Performed Performed but Not Available
## 2 2
## Positive
## 339
table(cdat$PR)
##
## Indeterminate Negative
## 3 147
## Not Performed Performed but Not Available
## 2 2
## Positive
## 291
table(cdat$HER2)
##
## Equivocal Negative Not Available Positive
## 5 370 4 65
table(cdat$Node)
##
## 0 1 2 3
## 221 146 54 23
table(cdat$Metastasis)
##
## 0 1
## 427 11
table(cdat$ER,cdat$PR)
##
## Indeterminate Negative Not Performed
## Indeterminate 0 1 0
## Negative 1 93 0
## Not Performed 0 0 2
## Performed but Not Available 0 0 0
## Positive 2 53 0
##
## Performed but Not Available Positive
## Indeterminate 0 1
## Negative 0 6
## Not Performed 0 0
## Performed but Not Available 2 0
## Positive 0 284
#cluster heatmap - biclustering
#cluster heatmap - biclustering
aa = grep("grey",colors())
bb = grep("green",colors())
cc = grep("red",colors())
gcol2 = colors()[c(aa[1:2],bb[1:25],cc[1:50])]
Without scaling
heatmap(gdat,col=gcol2,hclustfun=function(x)hclust(x,method="ward.D"))
With scaling
heatmap(scale(gdat),col=gcol2,hclustfun=function(x)hclust(x,method="ward.D"))
Cols=function(vec){cols=rainbow(length(unique(vec)))
return(cols[as.numeric(as.factor(vec))])}
heatmap(scale(gdat),col=gcol2,hclustfun=function(x)hclust(x,method="ward.D"),labRow=cdat$Subtype,RowSideColors=Cols(cdat$Subtype))
#Dimension Reduction
PCA
Cols=function(vec){cols=rainbow(length(unique(vec)))
return(cols[as.numeric(as.factor(vec))])}
sv = svd(scale(gdat,center=TRUE,scale=FALSE))
V = sv$v
Z = gdat%*%V
K = 3
pclabs = c("PC1","PC2","PC3","PC4")
par(mfrow=c(1,K))
for(i in 1:K){
j = i+1
plot(Z[,i],Z[,j],pch=16,xlab=pclabs[i],ylab=pclabs[j],col=Cols(cdat$Subtype))
}
legend(-45,0,pch=16,col=rainbow(5),levels(cdat$Subtype))
Pairs Plot
PC1<-as.matrix(Z[,1])
PC2<-as.matrix(Z[,2])
PC3<-as.matrix(Z[,3])
PC4<-as.matrix(Z[,4])
PC5<-as.matrix(Z[,5])
pc.df.cdat<-data.frame(Subtype = cdat$Subtype, PC1, PC2, PC3, PC4, PC5)
ggpairs(pc.df.cdat, mapping = aes(color = Subtype))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
MDS
Dmat = dist(gdat,method="maximum")
mdsres = cmdscale(Dmat,k=2)
plot(mdsres[,1],mdsres[,2],pch=16,col=Cols(cdat$Subtype), main = "Dimension Reduction MDS")
legend(-30,20,pch=16,col=rainbow(5),levels(cdat$Subtype))
ICA
require("fastICA")
## Loading required package: fastICA
K = 4
icafit = fastICA(gdat,n.comp=K)
kk = 3
pclabs = c("ICA1","ICA2","ICA3","ICA4")
par(mfrow=c(1,kk))
for(i in 1:kk){
j = i+1
plot(icafit$A[i,],icafit$A[j,],pch=16,xlab=pclabs[i],ylab=pclabs[j],col=Cols(cdat$Subtype))
}
legend(-1,2.8,pch=16,col=rainbow(5),levels(cdat$Subtype))
UMAP
gdat.umap = umap(gdat)
plot(gdat.umap$layout[,1], y =gdat.umap$layout[,2], type = "n", main = "UMAP", xlab = "UMAP1", ylab = "UMAP2")
text(gdat.umap$layout[,1], y =gdat.umap$layout[,2], type = "n", cdat$Subtype, col=Cols(cdat$Subtype), cex = .7 )
## Warning in text.default(gdat.umap$layout[, 1], y = gdat.umap$layout[, 2], :
## graphical parameter "type" is obsolete
#Clustering
K-means
K = 5
km = kmeans(gdat,centers=K,nstart=25)
table(km$cluster,cdat$Subtype)
##
## Basal-like HER2-enriched Luminal A Luminal B Normal-like
## 1 5 8 45 12 4
## 2 0 2 42 60 1
## 3 74 5 1 1 1
## 4 0 7 109 26 1
## 5 0 31 3 7 0
Plot Kmeans with labels
plot(Z[,1],Z[,2],col=km$cluster, main = "Plot Kmeans Clusters ", xlab = "PC1", ylab = "PC2")
text(Z[,1],Z[,2],cdat$Subtype,cex=.75,col=km$cluster)
cens = km$centers
points(cens%*%V[,1],cens%*%V[,2],col=1:K,pch=16,cex=3)
Hierarchical
#which linakge is the best?
#which distance metric is the best?
Dmat = dist(gdat)
com.hc = hclust(Dmat,method="ward.D")
plot(com.hc,labels=cdat$Subtype,cex=.5)
res.com = cutree(com.hc,5)
table(res.com,cdat$Subtype)
##
## res.com Basal-like HER2-enriched Luminal A Luminal B Normal-like
## 1 1 3 95 11 3
## 2 0 4 73 65 1
## 3 75 4 5 4 1
## 4 0 27 3 7 0
## 5 3 15 24 19 2
Consensus Clustering with Hierarchical
results = ConsensusClusterPlus(gdat,maxK=6,reps=500,pItem=0.8,pFeature=1,
clusterAlg="hc",distance="pearson",plot="png")
## end fraction
## clustered
## clustered
## clustered
## clustered
## clustered
Look at results for first 5 clusters
heatmap(results[[2]][["consensusMatrix"]][1:5,1:5])
Spectral Clustering
K = 5
s_gdat = specClust(gdat, centers=K, nn = 7, method = "symmetric", gmax=NULL)
Visualize
plot(Z[,1],Z[,2],col=s_gdat$cluster, main = "Visualize Spectral Clusters", xlab = "PC1", ylab = "PC2")
text(Z[,1],Z[,2],cdat$Subtype,cex=.75,col=s_gdat$cluster)
Genes significantly associated with ER or PR Status, etc
x = gdat[cdat$ER=="Positive" | cdat$ER=="Negative",]
y.er = cdat$ER[cdat$ER=="Positive" | cdat$ER=="Negative"]
y.label = rep(1, length(y.er))
y.label[y.er == "Positive"]=2
ps = NULL
for(i in 1:ncol(gdat)) ps = c(ps,
t.test(x[y.label==1,i],x[y.label==2,i])$p.value)
fdrs.bh = p.adjust(ps, method="BH")
cat("Number of Tests significant with alpha=0.1 using Bonferroni correction:",
sum(ps<0.1/length(y.label)), fill=TRUE)
## Number of Tests significant with alpha=0.1 using Bonferroni correction: 165
cat("Number of Tests with FDR below 0.1:",
sum(fdrs.bh<0.1), fill=TRUE)
## Number of Tests with FDR below 0.1: 259
plot(sort(ps,decreasing=FALSE),ylab="P-Values")
#BH procedure
abline(a=0, b=0.1/length(y.label),col="red")
#Bonferroni
abline(a=0.1/length(y.label), b=0,col="blue")
Graphical models - how are genes related?
#this takes a bit of time
lam = lambdaMax(gdat)*sqrt(log(ncol(gdat))/nrow(gdat))*0.01
net = XMRF(t(gdat),method="LPGM",lambda.path=lam,N=1,th=.0025)
## LPGM ::: Conducting sampling ... in progress: 100 %
## LPGM Completed.
plot(net)
## Plot optimal network
## [,1] [,2]
## [1,] 19.838264657 30.86689
## [2,] 8.295297756 30.98339
## [3,] -4.945591598 48.50826
## [4,] -0.072424102 38.72952
## [5,] 1.446459793 35.98709
## [6,] 3.194760332 39.58998
## [7,] 21.312783923 30.77726
## [8,] 8.373994622 35.18537
## [9,] 4.707045609 57.09021
## [10,] 5.158871953 33.55555
## [11,] -2.730617466 37.54865
## [12,] 5.678206097 56.22942
## [13,] 25.438684620 40.82255
## [14,] 8.291445725 55.49436
## [15,] 3.042894993 56.59438
## [16,] -1.263971134 50.76479
## [17,] 26.142476160 45.51535
## [18,] 12.425998424 31.38349
## [19,] 8.502924125 41.45538
## [20,] -0.200305414 33.08569
## [21,] 16.081662098 51.65025
## [22,] 21.917792208 51.52781
## [23,] 21.832471278 47.67047
## [24,] 26.599573683 41.47744
## [25,] 3.565047891 37.93286
## [26,] 4.339860804 31.75050
## [27,] 7.324784847 52.33422
## [28,] 20.228125422 34.05805
## [29,] 4.852517564 31.08728
## [30,] 3.292012143 31.96992
## [31,] -1.607282823 40.33197
## [32,] 11.220051034 35.32044
## [33,] -0.480258757 53.37978
## [34,] -1.903283851 37.86544
## [35,] 20.999739217 43.35891
## [36,] 23.182693792 33.01077
## [37,] 6.863407060 57.16108
## [38,] 9.367402156 44.12062
## [39,] 9.377171601 49.81023
## [40,] 18.671656093 56.29594
## [41,] 2.724928650 34.16188
## [42,] 7.612850869 41.95366
## [43,] 18.230046824 46.32492
## [44,] 0.660423176 38.95581
## [45,] 13.898747569 35.62331
## [46,] 13.841358863 56.75667
## [47,] 19.540657234 47.48261
## [48,] 1.343899473 32.35624
## [49,] 10.834075224 56.98855
## [50,] 22.109606361 53.52758
## [51,] -3.960265652 47.45734
## [52,] 1.569623673 56.45270
## [53,] 25.932057403 37.54140
## [54,] 1.544511485 37.29487
## [55,] 10.765269363 34.83927
## [56,] 20.225106604 48.75143
## [57,] 16.589417314 39.41186
## [58,] 11.102320523 36.93196
## [59,] 14.057238179 36.67035
## [60,] -3.445348623 51.58113
## [61,] 1.444924244 40.59768
## [62,] 10.489475429 32.71512
## [63,] 1.030163536 40.40300
## [64,] 9.296411878 34.64020
## [65,] 14.347754211 38.13310
## [66,] 22.678879007 40.97109
## [67,] 20.026789236 55.50454
## [68,] 13.034252033 55.40708
## [69,] 17.519550403 53.03269
## [70,] 19.483331214 32.10677
## [71,] 26.363384270 44.13467
## [72,] 13.419331348 36.90529
## [73,] 0.670777338 39.97927
## [74,] -3.683297526 49.08652
## [75,] 7.900145849 34.66698
## [76,] 2.802066170 38.97391
## [77,] -2.031177041 53.58478
## [78,] 4.737252968 36.27136
## [79,] 5.820065189 36.89361
## [80,] 0.904722613 31.75210
## [81,] -2.470066043 48.47391
## [82,] 5.853995445 50.76320
## [83,] 10.120748167 41.99061
## [84,] 13.868630933 58.21228
## [85,] 18.428870221 43.32143
## [86,] 1.007264220 52.35334
## [87,] 20.355687570 37.99708
## [88,] -0.196485013 42.40738
## [89,] 22.865013156 44.52788
## [90,] 4.295255982 40.56307
## [91,] 8.838605013 51.75406
## [92,] 0.798961138 54.12138
## [93,] 20.549367075 35.59592
## [94,] 11.364135124 30.94853
## [95,] 8.146185457 37.44643
## [96,] 5.765998187 39.80018
## [97,] 15.003234626 46.87180
## [98,] 7.534188860 50.54428
## [99,] 9.401946764 56.36881
## [100,] 10.515262743 35.83292
## [101,] 15.671676728 27.14843
## [102,] 11.398655079 55.64188
## [103,] 12.744248299 53.83570
## [104,] 1.792001410 38.67687
## [105,] 11.792116145 31.51406
## [106,] 7.784714063 31.11862
## [107,] 16.152423957 54.74328
## [108,] 21.971439597 32.44607
## [109,] 7.150080834 36.21841
## [110,] 2.463621634 51.29772
## [111,] 17.347555735 55.97895
## [112,] -1.971338480 36.86097
## [113,] 12.549164880 39.63456
## [114,] 11.156177403 58.63994
## [115,] 19.793871747 45.91049
## [116,] 14.290558433 37.73597
## [117,] 5.197463501 36.48672
## [118,] 24.671787425 50.02539
## [119,] -1.024601796 41.11683
## [120,] 1.546314270 39.72496
## [121,] 15.240984587 57.90895
## [122,] 17.516320397 50.84673
## [123,] 13.262042582 39.62001
## [124,] 3.558381286 36.67464
## [125,] 3.088762484 33.63772
## [126,] 13.751097631 51.41151
## [127,] 9.807747524 54.80725
## [128,] 4.658885303 44.19630
## [129,] 11.029461087 33.43757
## [130,] 4.504082787 33.56112
## [131,] 3.689740068 36.99346
## [132,] 4.016149566 42.37604
## [133,] 24.062312752 40.30369
## [134,] 7.282962721 58.66570
## [135,] 5.542570977 41.10222
## [136,] 5.164241054 31.68914
## [137,] 21.361638072 44.86243
## [138,] 10.443698628 52.73173
## [139,] 4.617545296 39.94560
## [140,] 1.805002569 41.21817
## [141,] -0.625038606 43.38235
## [142,] 5.922463717 53.17410
## [143,] -0.795226171 33.36801
## [144,] 9.505580830 36.05210
## [145,] 5.947356835 41.23146
## [146,] 6.741400740 43.88794
## [147,] 6.232160576 38.06550
## [148,] 18.174144724 29.47499
## [149,] 18.984864842 29.75727
## [150,] 8.194926084 57.56744
## [151,] 2.737755542 37.92951
## [152,] 3.566641382 40.61698
## [153,] 13.421756629 35.27450
## [154,] 12.961315673 38.91233
## [155,] 2.201081793 34.97813
## [156,] 14.300970591 34.29035
## [157,] 14.522505545 35.91764
## [158,] 14.194409995 35.07020
## [159,] 5.794803148 33.80872
## [160,] 23.048211645 39.16700
## [161,] 13.221070178 36.27409
## [162,] 2.173762497 53.82266
## [163,] 7.131777364 33.46496
## [164,] 22.602932716 31.35726
## [165,] -1.242489866 42.85625
## [166,] 24.184834115 38.28010
## [167,] 11.058047615 30.40826
## [168,] 16.914687337 28.00923
## [169,] 0.595773901 36.44548
## [170,] 8.727065121 36.19884
## [171,] 3.503497471 31.44125
## [172,] 23.008958362 52.38288
## [173,] 3.502070291 36.37646
## [174,] 1.088052448 38.82101
## [175,] 13.852374093 36.18074
## [176,] 12.271326498 35.70071
## [177,] 11.831185597 33.42746
## [178,] 6.627087544 34.57717
## [179,] -5.305666592 46.59879
## [180,] 13.528831284 47.98797
## [181,] 24.605389472 46.43136
## [182,] 17.020539887 57.39635
## [183,] 9.426264208 31.30048
## [184,] 1.355908408 35.27176
## [185,] 12.231785935 52.32198
## [186,] 9.420958674 38.06781
## [187,] 8.943182968 53.49179
## [188,] 26.394092588 39.62829
## [189,] 14.829038658 55.71897
## [190,] 6.832353517 55.45470
## [191,] -2.358970474 39.54544
## [192,] -1.666506659 34.04969
## [193,] -2.275242197 52.24315
## [194,] 0.339304584 50.77731
## [195,] 8.551931395 35.66876
## [196,] 21.490377333 39.61240
## [197,] 25.301002346 48.27629
## [198,] 4.008223791 50.72216
## [199,] 25.306556265 38.95048
## [200,] 8.860426188 30.77092
## [201,] -0.281277037 35.18854
## [202,] 0.539689078 35.22785
## [203,] 14.283730116 52.82004
## [204,] 16.446448589 46.19847
## [205,] 18.509402870 51.84716
## [206,] -1.287485890 36.21155
## [207,] 0.532260756 40.77345
## [208,] 8.780027709 31.58064
## [209,] 2.132271940 46.57556
## [210,] 2.905911066 44.07795
## [211,] -0.117219586 32.40754
## [212,] 9.109844984 37.12172
## [213,] 4.884630741 34.96444
## [214,] 3.990345452 35.60359
## [215,] 6.306044925 39.00490
## [216,] 20.121503415 51.39196
## [217,] 21.567968334 36.46362
## [218,] 9.346722849 33.57526
## [219,] 8.480829085 39.30396
## [220,] 1.805003731 37.88845
## [221,] 8.942536751 40.12915
## [222,] 12.327058669 57.07843
## [223,] 0.530713842 37.01683
## [224,] -0.202187719 37.73361
## [225,] 5.881059705 39.26675
## [226,] -0.434982039 36.91457
## [227,] 20.754545431 32.76185
## [228,] -2.549200370 50.29915
## [229,] 7.427969542 37.70748
## [230,] 13.794024250 39.06348
## [231,] 26.163350380 42.71873
## [232,] 20.882350744 52.73482
## [233,] 6.298537381 42.66674
## [234,] 3.075718331 35.74845
## [235,] 9.117956833 40.45187
## [236,] 22.289352099 42.49072
## [237,] 9.378959486 59.03277
## [238,] 0.461578119 32.35721
## [239,] 11.962596039 50.65777
## [240,] 1.629106732 36.40548
## [241,] 15.930469830 56.58868
## [242,] 17.370782612 44.81138
## [243,] 9.741241772 31.83634
## [244,] 15.422347606 32.64791
## [245,] 2.992743180 52.61110
## [246,] 13.324546472 49.61528
## [247,] 0.861906584 37.36720
## [248,] 3.746359467 35.07027
## [249,] 4.883874289 52.09552
## [250,] 3.829502191 33.80065
## [251,] 23.706780285 51.17331
## [252,] 10.578815917 34.12663
## [253,] 9.855902993 32.74570
## [254,] 10.398643425 30.19314
## [255,] 24.095915888 48.75595
## [256,] 25.840009426 46.96213
## [257,] -0.672700369 52.08543
## [258,] 5.450968511 33.94353
## [259,] 21.897124164 37.98763
## [260,] 9.961617224 35.69469
## [261,] 17.661452824 54.53311
## [262,] 14.504123072 36.45517
## [263,] 23.090094348 46.01576
## [264,] 3.196178189 42.98282
## [265,] 5.678252169 31.34060
## [266,] 23.421533470 47.53541
## [267,] -0.742544464 39.61887
## [268,] 21.102362475 50.20402
## [269,] -3.429137838 41.17751
## [270,] 13.320598350 34.63415
## [271,] 7.405313263 37.45012
## [272,] 2.705068806 39.73527
## [273,] 20.442421975 29.44126
## [274,] 15.968419945 53.14040
## [275,] 11.230467717 54.14501
## [276,] 3.942458003 55.26329
## [277,] -0.116401910 39.28005
## [278,] 2.783388445 36.34478
## [279,] 0.482561786 55.56206
## [280,] 19.563960474 42.13002
## [281,] 2.520830116 55.20942
## [282,] 2.590745759 35.24679
## [283,] -0.576073795 40.37840
## [284,] 9.797458107 40.55278
## [285,] 7.948746557 41.28381
## [286,] 7.765386982 40.88162
## [287,] 5.397004706 40.58520
## [288,] 4.986185106 33.19470
## [289,] 0.760662853 38.42472
## [290,] 24.655935356 36.79929
## [291,] 24.571517037 44.93621
## [292,] 3.076881019 41.44441
## [293,] 24.923961817 43.50148
## [294,] 2.505749609 36.85567
## [295,] 4.674507267 46.79393
## [296,] 21.028445051 41.44502
## [297,] 24.247606757 33.72549
## [298,] 1.542356920 50.18198
## [299,] 7.081825746 35.44002
## [300,] 17.159343282 47.79718
## [301,] 21.936724774 34.23455
## [302,] 9.661174753 57.83958
## [303,] 19.103697846 54.49899
## [304,] 14.416119074 54.26214
## [305,] -0.909443890 49.33533
## [306,] 18.400181200 28.08469
## [307,] 15.259638461 48.61540
## [308,] 5.279185558 35.66529
## [309,] 7.389919538 32.94221
## [310,] 19.362738935 53.09857
## [311,] 23.555292973 43.15174
## [312,] 23.949805001 35.01009
## [313,] 15.980224429 33.49990
## [314,] 19.268796328 50.07388
## [315,] 2.039220549 41.85566
## [316,] -0.786454720 54.79682
## [317,] 19.534304324 44.40725
## [318,] 16.622968016 49.59964
## [319,] 22.936077684 50.00550
## [320,] 11.523671387 49.00336
## [321,] -4.301039961 50.28031
## [322,] 21.355132041 46.40529
## [323,] 4.092743574 53.56517
## [324,] 3.354968173 38.87101
## [325,] 3.629821483 57.72785
## [326,] 8.759446640 34.10995
## [327,] 12.830559872 34.96710
## [328,] 13.550685049 38.41369
## [329,] -0.005711561 43.81178
## [330,] 10.431439722 51.06879
## [331,] 5.470333110 38.35883
## [332,] 7.756134211 38.58301
## [333,] 2.228544534 35.61018
## [334,] 25.191868878 35.58392
## [335,] 7.377755388 54.05618
## [336,] 20.873025165 54.34941
## [337,] 4.107299070 31.11167
## [338,] 24.403930085 41.89686
## [339,] 14.995256574 50.43140
## [340,] 12.514656907 58.61875
## [341,] 5.681890937 58.22659
## [342,] 22.714005711 35.44925
## [343,] 23.138558363 37.01331
## [344,] 20.030065198 40.11837
## [345,] 1.231276141 43.27957
## [346,] 4.466884527 38.17857
## [347,] 7.450070180 32.09569
## [348,] 1.960576018 39.02725
## [349,] 7.389100441 43.11813
## [350,] 18.299190485 48.85301
## [351,] 22.034390427 49.00946
## [352,] 9.916337207 38.30725
## [353,] 5.277049733 54.67881
## estimation using spacejam
net2 = SJ(gdat)
net2 = net2$G
sjnet <- 1*(net2)[,,1]
sjnet = graph.adjacency(sjnet, mode="undirected")
plot(sjnet, vertex.size=0.1, vertex.label=NA)